%To write the source current files 
% for cable computation
%latest date 29.9.5
MyInPath = 'c:\Manoj\ocean\Velocities\monthly\';
MyOutPath = 'c:\Manoj\ocean\Source\180x360\monthly\';

if isdir(MyInPath) == 0,
    display('In path does not exist');
    return;
end;
if isdir(MyOutPath) == 0,
    mkdir('c:\Manoj\ocean\','Source\180x360\monthly\');
end;

S = dir(MyInPath);

nfiles = length(S);
ndatafile=1;
for i = 1:nfiles,
    if length(S(i).name) == 8,
        FileName(ndatafile,:) = S(i).name;
        ndatafile = ndatafile+1;
    end;
end;
fp = fopen('projectfiles.txt','wt');
%for i = 1:ndatafile-1,
for i = 1:1,
    disp([MyInPath FileName(i,:)]);
 fid = fopen([MyInPath FileName(i,:)],'r','b');
 IC = fread(fid,[180,360],'float');% integrated conductivity - Not applicable in your case
 ICU = fread(fid,[180,360],'float');% integrated conductivity times IU- Not applicable in your case
 ICV= fread(fid,[180,360],'float');%integrated conductivity times IV- Not applicable in your case
 IU = fread(fid,[180,360],'float');%integrated velocity Zonal U, East 
 IV = fread(fid,[180,360],'float');%integrated velocity Meridional V, North 
 Bz = fread(fid,[180,360],'float');% Bz from a main field model (SHD upto 16 only required)
 fclose(fid);
 
CU = -IV.*Bz.*1e-9.*3.2; %J^ext , zonal
CV = -IU.*Bz.*1e-9.*3.2;%J^ext , meridional
CU = circshift(CU,[0,180]); % change the longitude from -180 180 to 0-360
CV = circshift(CV,[0,180]);

fprintf(fp,'%s\n',['    4000.           ' FileName(i,1:end-3) ' ' FileName(i,1:end-3) '.source' ' 2']); %open a .source file for writing
fid = fopen([MyOutPath FileName(i,1:end-3) '.source'],'wt');
fid1 = fopen('c:\manoj\ocean\x3dg\header.txt','rt');
%the following format for alexei's x3dg is tested ok 20 oct 05 - pl don't
%mess with them :-)
           for ii = 1:34,
                s=fgetl(fid1);
                fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:180,
    for jn = 1:360,
       fprintf(fid,'%12.4e\n',CU(in,jn)); %write real part of the J^ext zonal values 
   end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
for in = 1:180,
    for jn = 1:360,
        fprintf(fid,'%12.4e\n',0); %write imaginary part of the J^ext zonal values - in case of long period this is practically zero
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
 for in = 1:180,
    for jn = 1:360,       
        fprintf(fid,'%12.4e\n',CV(in,jn)); %write real part of the J^ext meridional values 
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:17,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);
  
 for in = 1:180,
    for jn = 1:360,       
        fprintf(fid,'%12.4e\n',0);%write imaginary part of the J^ext meridional values - in case of long period this is practically zero
    end;
end;
            fprintf(fid,'\n');
            for ii = 1:37,
            s=fgetl(fid1);
            fprintf(fid,'%s\n',s);
            end;
            s=fgetl(fid1);

           fclose(fid);
           fclose(fid1);
end;
fclose all;